rand('seed',1);
% spnet.m: Spiking network with axonal conduction delays and STDP
% Created by Eugene M.Izhikevich.                February 3, 2004
% Modified to allow arbitrary delay distributions.  April 16,2008
M=50;                 % number of synapses per neuron
D=5;                  % maximal conduction delay 
% excitatory neurons   % inhibitory neurons      % total number 
Ne=80;                Ni=20;                   N=Ne+Ni;
a=[0.02*ones(Ne,1);    0.1*ones(Ni,1)];
d=[   8*ones(Ne,1);    2*ones(Ni,1)];
sm=10;                 % maximal synaptic strength

% Take special care not to have multiple connections between neurons
s = cell(D,1);
for i = 1:D
    s{i} = sparse(N,N);
end

excitatory_type = 1;
inhibitory_type = 2;
syndata(excitatory_type).min = 0;
syndata(excitatory_type).lrate = 1;
syndata(inhibitory_type).min = -10;
syndata(inhibitory_type).lrate = 0;

for n=1:Ne
    p=randperm(N);
    k=1;step= floor(M/D);
    for i = 1:D
        s{i}(n, p(k:(k+step-1)))=8*ones(1,step);
        
        syntype{i}(n, p(k:(k+step-1)))= excitatory_type*ones(1,step);
        k=k+step;
        % delays?
    end;
  
end



for n=Ne+1:N
    p=randperm(Ne);

    s{1}(n,p(1:min(M,floor(Ne/2)))) = -5*ones(1,min(floor(Ne/2),M));
    syntype{1}(n,p(1:min(M,floor(Ne/2)))) = inhibitory_type*ones(1,min(floor(Ne/2),M));

    %post(i,1:min(Ne,M))=p(1:min(Ne,M));
    %delays{i,1}=1:min(Ne,M);                    % all inh delays are 1 ms.
end;

sd = cell(D,1);
for i = 1:D
    sd{i} = spones(s{i})*.01;
end

%s=[6*ones(Ne,M);-5*ones(Ni,M)];         % synaptic weights

ltp_decay = exp(-1/20)
ltpcurve = .1*cumprod([ 1 ltp_decay*ones(1,49)]);
ltd_decay = exp(-1/20)
ltdcurve = -.12*cumprod([ 1 ltd_decay*ones(1,49)]);
stdpcurve = [ fliplr(ltpcurve)  ltdcurve 0];
  

v = -65*ones(N,1);                      % initial values
u = 0.2.*v;                             % initial values
firings=[];                         % spike timings
firehist=cell(N,1);

I=zeros(N,D);
for sec=1:60*60*24                      % simulation of 1 day
  for t=1:1000                          % simulation of 1 sec
    
    randindex = ceil((N-3)*rand);
    randindex = randindex:randindex+2;
    I(randindex,1)=I(randindex,1)+20;                 % random thalamic input 

    fired = find(v>=30);                % indices of fired neurons
    spfired = sparse(v'>=30);
    v(fired)=-65;  
    u(fired)=u(fired)+d(fired);
    
    % deliver all the synapses to post targets with correct delays.
    % I is a shift register for each neuron.
    for i = 1:D
        I(:,i) = I(:,i) + (spfired*s{i})';
    end
        
    v=v+0.5*((0.04*v+5).*v+140-u+I(:,1));    % for numerical 
    v=v+0.5*((0.04*v+5).*v+140-u+I(:,1));    % stability time 
    u=u+a.*(0.2*v-u);                   % step is 0.5 ms

    
    I= [I(:,2:D) zeros(N,1)];  % rotate.     
    
    % save the spike data.
    firings=[firings;t*ones(length(fired),1),fired];

    for i =1:length(fired)
        firehist{fired(i)}(end+1)=t;
    end
    %imagesc(I);
  
    %drawnow;
  
  end;
  sec
  
  plot(firings(:,1),firings(:,2),'.');
  axis([0 1000 0 N]); drawnow;
  
  
  % do some learning based on spike times.
  % sd in Eugene's spnet did not get used until after 1 second, so this is
  % equivalent.  Just makes for a cleaner implementation.
  for delay = 1:D
    % get all the nonzero elements from synaptic matrix.
    [pre_i post_i] = find( s{delay} );
    
    % Adjust these weights.
    for j = 1:length(pre_i)
        type = syntype{delay}(pre_i(j),post_i(j));
        sd{delay}(pre_i(j),post_i(j)) = sd{delay}(pre_i(j),post_i(j)) + stdp(firehist{pre_i(j)},firehist{post_i(j)}, delay, stdpcurve,syndata(type).lrate);
        s{delay}(pre_i(j),post_i(j)) = max( syndata(type).min ,min(sm, syndata(type).lrate*.01 + s{delay}(pre_i(j),post_i(j)) + sd{delay}(pre_i(j),post_i(j)))); 
    end
    
    sd{delay}=0.9*sd{delay};
    
  end
  % keep spikes near the seconds boundary for future STDP calculations.
  for i = 1:N
    firehist{i}=firehist{i}( find(firehist{i}> 950) )-1000;
  end
  
  firings=[];
  
end;


